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Abstract 

Using molecular dynamics(MD) simulations, we study thermal transport in graphene nanorib- 
bons (GNR) subjected to uniform uniaxial and nonuniform strain fields. We predict significant 
thermal rectification (over 70%) in a rectangular armchair GNR by applying a transverse force 
asymmetrically. The heat flux is larger from the less stressed region to the more stressed region. 
Furthermore, we develop a theoretical framework based on the non-equilibrium thermodynamics to 
discuss when thermal rectification under a stress gradient can occur. We conclude with a discussion 
of details relevant to experiment. 
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I. INTRODUCTION 



Operation of nanoscale thermal devices, mainly relies on the tunability of phonon trans- 
port by external means. Thermal devices, such as thermal rectifiers^, thermal transistors^ 
and thermal memories^ are a new class of devices, whose operation is driven by the tem- 
perature gradients. These devices will have useful applications not only in thermal circuits 4 
but also in nanoscale thermal management and thermo-electric applications. 

A thermal rectifier, in which the thermal current is larger in one direction than in 
the opposite, is one of the most fundamental thermal devices to be realized experimen- 
tally in nanoscale. Thermal rectification has been demonstrated using molecular dynamics 
simulations (MD) in asymmetric systems and is discussed as resulting from an interplay 
between structural asymmetry and lattice anharmonicityi^. Recently, this phenomenon 
has been observed experimentally in asymmetrically mass-loaded carbon and boron nitride 
nanotubes 7 . Moreover, thermal rectification has been predicted in triangular and trape- 
zoidal shaped graphene nanoribbons(GNRs)^. In this letter, we introduce a method based 
on strain engineering as an experimentally feasible approach to realize significant (> 70%) 
thermal rectification in GNRs. The potentially real-time tunability of the thermal rectifica- 
tion is a distinctive advantage of this approach. 

Graphene is a promising material for nanoscale applications due to its exceptional 
electronic^, thermal and mechanical properties. It has been shown recently that the elec- 
tronic properties of graphene can be tuned significantly by engineering strain. For instance, 
strain could produce a pseudo-magnetic fielcr^ - — which affects the electronic properties of 
graphene. Furthermore, strain could induce semiconducting properties on metallic GNRs by 
opening a gap in the electronic band structure^ 4 -. However, strain induced tuning of thermal 
transport properties of graphene has been less studied. The effect of uniform uniaxial strains 
on the thermal conductivity (n) of graphene has been studied recently by MD simulations 
and the reported effect is substantial in tuning 

When increasing the uniaxial tensile strain (arrows in figJT^i), we observe broadening the 
gap in the phonon dispersion of armchair GNR (AGNR) whereas the effect on the phonon 
dispersion of zigzag GNR (ZGNR) is mild. The calculated k of AGNR also reduce consid- 
erably when increasing the tensile strain compared to that of ZGNR. We predict thermal 
rectification on an asymmetrically stressed rectangular GNR. In our system, asymmetric 
stress is achieved by exerting a transverse force (F y ) on the top and bottom edge atoms of 
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FIG. 1. Schematic of the 8-ZGNR (a) and 11-AGNR (b). The triangles represent the thermostated 
atoms and the full squares are the fixed atoms. The horizontal arrows in (a) shows the direction 
of uniaxial tensile strain. In (b) the vertical arrows show the applied constant force F y on the top 
and bottom atoms of the right half of the AGNR. The horizontal arrows in (b) represent the stress 
developed in length (x) direction (a R and cx£) near the heat baths. 

part (for example, the right half) of the AGNR as illustrated in fig{]Tj. We observe that the 
thermal current is larger from less stressed region to the more stressed region (left to right) 
than that in the opposite direction (right to left). Furthermore, we develop a theoretical 
frame work based on the non-equilibrium thermodynamics, to model thermal rectification 
in the presence of a stress gradient. 

II. SIMULATION PROCEDURE 

Figured] shows a schematic of n-ZGNR (Zigzag GNR) and n-AGNR structures simulated. 
The carbon-carbon bond length is 1.42 A in the absence of strains. The n refer to the 
number of carbon dimers ( a dimer is shown by the oval) in width 16 . The first and the 
last columns of atoms are fixed and the adjacent three columns of atoms of the n-ZGNR 
and four columns of atoms of the n-AGNR are coupled to the Nose -Hoover thermostats. 
In the MD simulation, a second generation Brenner potential is employed to describe the 
carbon-carbon interactions. The equations of motion are integrated using the third-order 
predictor-corrector method. The time step is 0.5 fs and the total simulation time is 5 ns 
(10 7 time steps). The temperatures of the left and right heat baths (HBs) are set to Tl 
and Tr respectively. The temperature difference (Tl — Tr) = 2aT, where a determines the 
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FIG. 2. Normalized thermal conductivity (k/kq) as a function of uniaxial strain (I — along 
the transport direction (x) calculated at 300 K. 

temperature bias and the average temperature T = (Tl + Tr)/2. The net heat flux (J) was 
calculated by the power delivered by the heat baths 18 . The thermal conductivity (k) of the 
system is calculated according to Fourier's law, J = nwhVT, where w and h are width and 
van der Waals thickness (h = 0.335 nm) of the GNR. 




III. THERMAL TRANSPORT UNDER UNIFORM STRAIN FIELDS 

Figure [2] shows the variation of k/kq of ZGNRs and AGNRs of two different widths 
as a function of strain ranging from —0.03 to 0.12. The k,q is the thermal conductivity of 
unstrained nanoribbon. The horizontal arrows in figfiji show the direction of uniaxial tensile 
strain on both nanoribbons. The strain is given by (I — Iq)/Iq, where I is the stretched or 
compressed length and Iq is the initial length of the sample. Strain is applied by slowly 
moving the fixed atoms on one side (right side) at a rate of ±1.0 x 10~ 6 A per time step 
at 300 K. The calculated values of the thermal conductivity (kq) of unstrained 8-ZGNR, 
12-ZGNR, 11-AGNR and 19-AGNR are 5360 W/mK, 5500 W/mK,3300 W/mK and 3600 
W/mK respectively. The variation of k for ZGNRs in the range of strain — 0.1 is about 30% 
from the unstrained values. The reduction of k at a strain of 0.1 is about 55% for 11-AGNR 
and 19-AGNR. Figure [3] shows the calculated phonon dispersions of 8-ZGNR and 11-AGNR 
along the transport direction(x) of unstrained nanoribbons and at a strain of 0.024. When 
increasing the uniaxial tensile strain, we do not observe a significant alteration of dispersion 
curves of the 8-ZGNR (figl3]a,b). However, we observe an increased of the phonon velocity 
(duj/dk) of some modes in the low (0-600 cm -1 ) and high (1250-1600 cm— 1) frequencies 
which could be the reason for the increment of the k in small strains. In 11-AGNR we 
observe broadening the narrow gap at 1350 cm -1 in the phonon dispersion. At a strain of 
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FIG. 3. Phonon dispersion of 8-ZGNR and 11-AGNR under uniaxial tensile strains. The horizontal 
axis is the normalized wave vector projected on to the transport direction(x). 



0.024, we observe a gap in the frequency range: 1250 — 1350 cm^ 1 (figJSJi). The gap further 
broaden and moves toward low frequencies when increasing the uniaxial strain. This could 
be the reasons for the relatively larger reduction of k of AGNR with increasing strain. 



IV. THERMAL TRANSPORT UNDER NONUNIFORM STRAIN FIELDS 

In many of the theoretical studies on thermal transport, the stress along the nanoribbon 
is considered to be uniform and uniaxial. We consider a situation that the stress at the left 
and right edges of the nanoribbon are o L and a R such that a L ^ a R . In the intermediate 
region it is assumed to vary smoothly. This kind of stress profile could possibly be realized 
experimentally using the differential thermal expansion of graphene grown on structured 
substrates^>2&2i or through externally applied strain 22 . We achieve such an asymmetric 
stress profile in our simulation by transversely stressing a half of the GNR. For example, we 
apply a constant force F y on the atoms of the top and bottom edges of the right half of the 
11- AGNR as depicted in figEEb. The tensile stress developed across the y direction of the 
right half is denoted by a v R . This lateral stress produces a tensile strain along the y direction 
and, hence a compressive strain along the x-direction due to the Poisson contraction. Since 
the atoms in the left and right edges are fixed, this compressive strain also results in tensile 
stress in the x-direction. The stress in the x-direction near the left and right HB is denoted 
by a x L and a x R . This is calculated from the time average of the forces on the fixed atoms at 
the left and right edges (excluding the corner atoms). Figure shows the variation of a x L 
and a R on F y calculated at T = 300 K. The a% and a R are tensile, and <r£ is smaller than 
the a R while both are found to be increasing with F y . 
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FIG. 4. (a) Variation of stresses {a R and <r£) on the 11-AGNR as a function of F y . These are 
the stresses in the center part of the nanoribbon near the right and left heat baths, (b) The ratio 
<x^/ct£ as a function of F y . The units of F y is nano Newton (nN). 
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FIG. 5. (a) Thermal current from left to right and right to left as a function of F y at T = 300 K 
and \a\ = 0.3, and (b) Variation of thermal current with the temperature bias (a) at F y = 4.8 nN. 



Thermal rectification 



In figEt) the thermal current from left to right (+a) and right to left (—a) are plotted 
against the applied force F y , when \a\ = 0.3 and T = 300 K. The thermal current from 
left to right {Jl-+r) shows very small decrement, where as the Jr^l decreases significantly 
when F y is increased. Thus, the Jl^r is considerably larger than Jr->l whenever F y > 0, 
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marking the existence of thermal rectification in this system. 

The tensile stress and strain can affect the thermal transport by distorting the lattice 
to both alter the characteristic vibrational frequencies and the degree of anharmonicity.— 
Anharmonicity is essential: no amount of geometric or parameter asymmetry can produce 
thermal rectification. Without anharmonicity, the system can be analyzed via normal modes, 
and time reversal symmetry will require that transmission amplitudes from left to right are 
the same as right to left. We attribute the observed thermal rectification to the strain 
induced asymmetry of the vibrational frequencies and the lattice anharmonicity. Such an 
asymmetry of the vibrational properties leads to a local variation of k. It can be shown that 
the thermal conductivity must be a function of both position and temperature to rectify the 
thermal current.— 

We can invoke the framework of non-equilibrium thermodynamics to describe this effect. 
The change of the entropy per unit volume of a solid due to the applied stress (r a ) can be 
expressed^, ds = ^ — A J2 a T adr] a , where du is the change of internal energy density, r] a is 
the strain and a 6 [xx, yy, xy\. From this we can deduce the rate of production of entropy^, 

s = ft (±) t ] + ±ft (r„) 4^ (1) 

where j\ and J^) are the z-th components of the energy and the strain currents, and 
repeated indices are summed over. In general the currents are a function of the intensive 
parameters (T and r) as well as the affinities^ (V|; and Vr). Thus the heat current, 
can be expanded in its most general form: 



jf) = L (Q) oL— + dj r a + L?.f dj-d k - + 

% l, J J rp ' rp J " ' lj,k J rp ^ rp ' 

dj-d k T a + ^ djT a d k T P (2) 

In solids, barring plastic deformation, there is no heat current in steady state solely due to 
Vrj. Thus, the coefficients Lfo a ' = L^^^ = 0. Moreover, in symmetric systems there is no 
thermal rectification, which implies L^®> = 0. If we further assume only gradients in the 
x-direction and stress in the y-direction the above equation can be reduced to, 

piQVvv) j 

jj^ = L^d x — H — d x —d x Tyy, (3) 
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FIG. 6. Thermal rectification factor (TR) at T = 300 K as a function of F y . The open circles 
are in the presence of edge disorder (EDO) and the open triangles are in the presence of vacancy 
defects (VD). 

which leads to different thermal currents when switching the sign of the stress gradient (Vr). 
For the case of F y = 4.8 nN, the calculated kinetic coefficients are LxQ = 1-7 x 10 8 WK/m 
and L& w) = 0.34 WK 2 /GPa. 

The distinctive property of this approach is the tunability of the thermal rectification by 
an asymmetrically applied force, F y . We show in figj6]the thermal rectification factor (TR) 
as a function of F y calculated at T = 300 K. The TR is defined as, 

TR = 2 \ Jl ^ r ~ Jr ^ l \ x 100%. (4) 
[Jl^r + Jr-^l) 

We observe an increment of the TR when increasing the F y . The increment of TR with F y 
is more prominent at higher biases (a = 0.3). As shown in the figj5b, the Jl->r increases 
almost linearly with the bias (a), whereas Jr^l increases nonlinearly giving rise to a larger 
TR. This behavior suggests the contribution of higher order terms (beyond the second order) 
to eqEJ The nonlinear transport is an essential property in realizing thermal rectification 
as in an electronic diode. In our system the nonlinear transport is prominent at larger F y s 
(> 0.8 nN), where we observe a clear bias dependence of TR. The maximum TR we observed 
is about 73% which occurred at a = 0.3 and F y = 5.6 nN. 

The right side of the AGNR is subjected to a biaxial stress whereas the left side only 
has an uniaxial stress. It is evident that the asymmetry of the on axis stresses (a R and 
cr£) is important in determining the TR. We achieve a significant TR whenever the ratio 
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cr^./o"£ > 2.0 (see figJU))- When we move the laterally stressed window towards the center, 
the TR reduced considerably in our simulation. This is a reasonable observation, because 
the asymmetry of the system is reduced. For instance, when we move the laterally stressed 
window eight columns towards the center (keeping the number of atoms that F y is applied 
constant) the TR decrease to 5% at a = 0.3 and F y = 4.8 nN. In this situation, the ratio of 
stresses reduces to 1.5 which is close to the ratio at F y = 0.8 nN where the observed TR is 
about 6%. When the length of the nanoribbon is increased, the strain gradient reduces and 
results a reduction of thermal rectification (implies from the eq]3j). However, by increasing 
the length of the AGNR and also the number of atoms that F y is applied, we could increase 
the strain gradient and observe a larger TR. 

Graphene commonly possesses edge disorder which significantly degrade the thermal 
properties. In figj6] the open circles shows the variation of TR in the presence of edge 
disorder. Since the simulated AGNRs are very small, we introduce only about a 4%( per- 
centage of number of edge atoms removed) edge disorder. We observe moderate reduction 
of TR, down to a value ~ 30%. The effect of edge disorder can be reduced by increasing 
the width of the nanoribbon. However, the presence of vacancy defects significantly reduces 
the TR. The open triangles in figj6] is for 0.6% vacancy defects distributed evenly in the 
nanoribbon. 

Experimentally, the lateral force on the right half of the AGNR can be applied by coupling 
to a substrate. Our simulation fixed the atoms to which the forces were applied. This 
situation also produces larger thermal rectification (over 100%) at higher biases and the 
direction of the maximum thermal current is same as before. This effect can be understood 
by the argument in ref.— . Consequently the observable net thermal rectification could be 
even higher due to the asymmetric coupling between the GNR and the substrate. 

Electronic transport of heat will occur in parallel to the phonon conduction, but is not 
so large that it dwarfs the phonon channel considered here. In addition, AGNRs shows 
semiconducting characteristics and their energy gap can be tuned with the strain^. Thus, 
the electronic contribution to the thermal transport should not be crucial. Electron-phonon 
interactions could lead to processes that undermine the thermal rectification. We believe 
that their contribution to the thermal rectification is also minimal since the nanoribbon is 
semiconducting. In addition, the long electronic coherence length in graphene indicates that 
electron-phonon interactions should not be significant. In experiments, these effects can be 
further minimized by electronically gating the sample. 
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Finally, we discuss the stability of the C-C bonds in graphene under the lateral forces. 
We do not observe any rupture of bonds within the applied range of F y (0 — 5.6 nN). In our 
simulation the estimated maximum force in the C-C bond direction is about 6.4 nN, which 
is correspond to a lateral stress of a y R = 90 Gpa. In a recent experiment, it has been found 
that the intrinsic strength of a single layer graphene is about 130 GPa assuming the van der 
Waals thickness of graphene^. When the F y is increased to 7.2 nN, we observe rupturing of 
some bonds near the fixed atoms in the right side of the AGNR. At this point, the maximum 
force in the C-C bond direction is about 8.0 nN. Thus, the forces required to realize thermal 
rectification by our method, are realistic and in an experimentally feasible range. 

V. CONCLUSION 

In conclusion, we study the thermal transport properties of strained GNRs using MD 
simulations. We demonstrate that the thermal rectification can be realized by engineering 
the stress on a rectangular AGNR. We have found that the heat transport is favorable from 
the less stressed region to the more stressed region. We also found that edge defects and 
vacancies only partially suppress this rectification. The major advantage of this approach 
is that the thermal rectification can be tuned from no rectification state to over 70% in 
real-time by applying a mechanical force. 
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